Growing point-to-set length scale correlates with growing relaxation times in model 

supercooled liquids 
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It has been demonstrated recently that supercooled liquids sharing simple structural features 
(e.g. pair distribution functions) may exhibit strikingly distinct dynamical behavior. Here we show 
that a more subtle structural feature correlates with relaxation times in three simulated systems 
that have nearly identical radial distribution functions but starkly different dynamical behavior. In 
particular, for the first time we determine the thermodynamic "point-to-set" length scale in several 
canonical model systems and demonstrate the quantitative connection between this length scale and 
the growth of relaxation times. Our results provide clues necessary for distinguishing competing 
theories of the glass transition. 
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PACS numbers: 61.43.Fs, 61.20.Lc, 64.70.Q-, 02.70.Ns 

The search for growing length scales accompanying the 
glass transition has been a major focus of the field for two 
decades Q. Great progress has been made in quantify- 
ing the behavior of dynamical length scales associated 
with emergent dynamical heterogeneity, such as £4, the 
length scale of the 4-point susceptibility . Less mature 
is our understanding of possible thermodynamic length 
scales that grow upon supercooling. In some theories 
of the glass transition, such as the random first-order 
theory (RFOT), dynamical arrest is connected to a par- 
ticular thermodynamically-based length scale associated 
with the depletion of independent particle configurations 
constrained by neighboring particles 0, HI • Within this 
viewpoint the putative structural length scale is distinct 
from £4 for modest supercooling Q • In other approaches 
this length scale may be identical to that associated with 
dynamical heterogeneity at all temperatures 0, Q . Fur- 
thermore, in these contrasting viewpoints the relation- 
ship of the various growing length scales with growing 
relaxation times are different. Clearly, elevating our level 
of knowledge of growing thermodynamic length scales as- 
sociated with glass formation to that of dynamical ones 
is of paramount importance in the continued quest for 
a deeper understanding of the behavior of supercooled 
liquids and glasses. 

A precise definition of one type of non-trivial struc- 
tural length scale akin to that envisioned in theories like 
RFOT was put forward by Bouchaud and Biroli in 2004 
Q . Their work suggests a procedure for the extraction of 
this length scale in computer simulations. This length is 
a measure of the distance scale over which particles are 
self-consistently pinned by other particles in their vicin- 
ity. The static correlations embodied by this length scale 
are known as "point-to-set" (PTS) correlations. In what 
follows, we will use £pts to denote the PTS length scale. 
There has been a recent flurry of activity in the extrac- 
tion of variants of £pts in a host of model glass forming 
systems 0, @-[3] ■ In this work we use this technique to 
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FIG. 1. (a) The radial distribution functions for A-type 
particles at T = 0.6 for the three models studied as an illus- 
tration of the structural similarity between the models. Other 
temperatures and pairs are shown in the supplemental infor- 
mation (SI), (b) The self-intermediate scattering functions 
for yl-type particles at k = 7.25 for NVE molecular dynamics 
simulations of the bulk systems. The alpha relaxation, r Q , 
is defined as the 1/e time of F s (k,t), illustrated by the hor- 
izontal dotted line. Curves shown for all models are, from 
L-R, T = {2.0,1.0,0.8,0.6,0.55} and additionally for WCA 
T = {0.5,0.45,0.4}. Note the LJ and IPL match at all tem- 
peratures while the WCA matches only at T=2.0 and decays 
significantly faster at lower temperatures. These temperature 
represent the regime computationally accessible to our cavity 
simulations. Parameters for the models can be found in the 
SI. 



address the crucial question of correlation between £prs 
and dynamics in systems with identical simple structural 
features (e.g. radial distribution functions, g(r)) but no- 
ticeably distinct relaxation time scales. 

We consider three closely related systems, chosen 
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specifically to address this problem. The first two are 
the standard Kob-Andersen binary Lennard- Jones mix- 
ture (LJ) [l5[ and its Weeks-Chandler- Andersen trunca- 
tion (WCA) [lfjj], which were previously found to have 
significantly different dynamical behavior at supercooled 
temperatures jl7j . despite having nearly identical two- 
body static correlations at all temperatures (e.g. mea- 
sured by g(r)). The third is a system characterized by a 
repulsive inverse power law potential (IPL) and was con- 
structed based on the composition and parameters of the 
LJ system to reproduce both two-body structural fea- 
tures and the dynamics of the LJ model [l8|. The sim- 
ilarities in pair correlations are illustrated in Fig. QJa), 
and the structural relaxation times for the temperatures 
in this study are shown in Fig. [ljb). While it should 
be noted that minor differences exist in g(r) between the 
three systems [HI (see supplemental information (SI)), 
they are far too small to account for the significantly 
weaker temperature dependence of relaxation times in 
the WCA system. In particular, neither mode-coupling 
theory nor activation-based theories that rely solely on 
the radial distribution function can account for this dis- 
tinction [ijj]. These three systems will allow us to study 
whether £pts is sensitive to structural differences in the 
three models not appearing at the level of g{r), and if so 
whether £pts can distinguish the LJ and IPL from the 
WCA. 

To measure £pts , we follow the protocol of Ref. [Io|- All 
simulations were perf ormed by the Monte Carlo (MC) 
technique of Ref. l2Cl Bulk equilibrium configurations 
are generated at a desired temperature, and then cavi- 
ties are constructed by freezing the particles outside a 
sphere of radius R. The center of the cavity is par- 
titioned into N cubes of side length I. The overlap is 

then defined as q(R,t) = (^N)- 1 J2?=i («i(*oK(*o + t)) 
where brackets denote both a thermodynamic average 
and an average over independent cavities, and rii(t) is a 
binary digit specifying whether a particle is in box i at 
"time" t (time here standing for any measure of simula- 
tion progress). As the cavity evolves in simulation, this 
quantity decays to a plateau which is independent of time 
and denoted q(R) = lim^oo q(R, t). Overlaps generated 
in this fashion will be referred to as "standard overlaps" . 
With the preceding definition, two independent config- 
urations will have an overlap q — pi 3 which is also the 
value of q(R —> oo). We will henceforth subtract off this 
bulk overlap from q, and denote the resulting value q. 
Further details can be found in the SI. 

A major complication of the algorithm sketched above 
is that the plateau value q{R) will over-estimate the true 
thermodynamic value of the overlap if the particle con- 
figuration inside the cavity breaks ergodicity or if the 
confinement simply induces relaxation on a time scale 
beyond that accessible to our simulations. We note that 
for all three systems studied here, relaxation times in- 



crease dramatically as cavity radius is decreased. The 
technique of particle swapping [2(| ameliorates this prob- 
lem in some systems, but is not effective in the systems 
studied here, as swap moves that exchange particles of 
different species are almost never accepted at supercooled 
temperatures. To test for convergence to the true (ther- 
modynamic) value of q(R), Cavagna and coworkers have 
proposed a technique based on the insertion of a random 
configuration of the same particles in the same cavity 
2l|. Initially such a configuration will have, on average, 
zero overlap (q — 0) with the pre-randomized configura- 
tion, but the structure of the boundary will induce finite 
overlap at long times. If this value measured with re- 
spect to the initial configuration yields the same value of 
q(R) as that extracted from the direct decay of q(R,t), 
one can be confident that the true thermodynamic value 
of the overlap has been obtained [2l[. This quantity is 
difficult to converge in our systems, hence we have imple- 
mented an approach which we call "particle size anneal- 
ing" (PSA). In this method the particles inside the cavity 
are instantaneously reduced in size such that their posi- 
tions quickly randomize, and are then evolved with the 
constraint of the cavity in place while their diameters are 
slowly tuned back to their original size. We have found 
this method, which is similar in spirit to algorithms used 
to g enerate randomly jammed packings of hard spheres 
[22| , is more efficient and reliable in the generation of con- 
verged overlap values than standard Monte Carlo sam- 
pling. In what follows, we use PSA both as a check of the 
convergence of q(R) to its thermodynamic value and as a 
means of generating estimates of q(R) for small R values. 
Specifically, for radii where standard and PSA overlap 
values match within error bars, the standard overlap is 
taken as the thermodynamic overlap. For smaller radii, 
the value from PSA yields a lower bound to the converged 
thermodynamic value of q(R). Further, we expect PSA 
to yield estimates of q(R) that are extremely close to the 
desired thermodynamic values. This expectation arises 
from a comparison with other sampling techniques (e.g. 
replica exchange) and direct Monte Carlo sampling. A 
discussion of these comparisons will be made in a future 
publication. Details of our approach may be found in the 
SI. Example results for a single cavity size can be seen in 
Fig. [Ha). 

We now turn to a discussion of the extraction of £pts 
from the spatial decay of the converged overlap function 
q(R). Overlaps were fit to a generalized compressed ex- 
ponential of the form q(R) = ^4exp(— ( g^~ J )■ Though 
a previous study used a pure compressed exponential (i.e. 
a = 0)[l0||, we choose a = 1, physically motivated by the 
fact that cavities with R ~ 1 should on average contain 
a single particle and that the overlap properties at this 
cavity size should not be sensitive to growing amorphous 
order. Furthermore we do not expect the same com- 
pressed exponential form to extend to cavities contain- 
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FIG. 2. (a) Overlap as a function of MC sweeps is shown 
for cavity size 72 = 3.1 at two different temperatures. Filled 
symbols show the overlap from MC dynamics and solid lines 
show stretched exponential fits to this data. Open symbols 
show the overlap from PSA, and dashed and dotted lines show 
horizontal fits to the long time plateau. Note that the overlap 
from PSA and MC dynamics do not meet at the lower tem- 
perature. In PSA, the particles have reduced diameters and 
may sample many more configurations than when they are 
full-sized, thus the PSA data go below the bulk value at short 
times. Error bars are from a bootstrap analysis and for the 
higher temperature are substantially smaller than the symbols 
shown, (b), (c) Overlap as a function of cavity size at two 
different temperatures. Here, closed symbols show standard 
overlaps, open symbols overlaps from PSA. Lines through the 
data are compressed exponential fits with the form discussed 
in the text. The LJ system is represented by solid lines, the 
WCA dashed, and the IPL dash-dotted. The insets show the 
same LJ data from the main figure with the same axis limits. 
Overlaid are data from the WCA system at which the length 
is most similar — T = 0.6 in (b) and T — 0.4 in (c). Boot- 
strap errors are, at largest, the size of the symbols shown. 



ing on average fewer than one particle. Fixing an a > 
allows us to perform a two-parameter fit with fixed A, 
leading to values of £pts with much smaller statistical 
variance. We find that for a — 1, A — 0.5 gives good 
fits to the data for all three systems at all temperatures 
studied; some example fits can be seen in Fig. [2jb) and 
HKc). An extended discussion of our fitting choices and 
methodology can be found in the SI, as well as a table of 
fit parameters extracted from the data. 

Fig. [3Ja) illustrates our first main result, namely the 



growth of the absolute thermodynamic length scale £pts 
for all three systems discussed above. Several notable fea- 
tures deserve mention. First, it is clear that the length 
scale £,pts grows unambiguously as temperature is low- 
ered. This is fully consistent with other recent studies 
that demonstrate growth of £pts in a variety of pinning 
geometries @,&[lJ|- Second, the distinction between the 
magnitude of £pts in the WCA system compared to the 
LJ and IPL systems at the same absolute temperature is 
stark. Despite the fact that the pair distribution func- 
tions of all three systems are nearly identical, the more 
subtle structural marker £pts can clearly distinguish the 
WCA system from the other two. The lengths of the IPL 
and LJ are found to be nearly identical at all tempera- 
tures and much larger than those found for the WCA. 

We now address the crucial question of correlation with 
relaxation times in a quantitative manner. Is £pts cor- 
related in a one-to-one manner with the alpha relaxation 
time r a extracted from the self-intermediate scattering 
functions of the systems under investigation? A key com- 
ponent of the answer to this question may be found in 
Fig. [3](b). While the absolute magnitude of £pts in the 
WCA system is clearly smaller than that of the other two 
systems at the same temperature, the lengths of all three 
systems collapse when temperature is scaled to the value 
where £pts ~ 1-4. These temperatures are quite similar 
to values obtained for "onset " temperatures obtained by 

In Fig. (3{c) we 
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independent means in earlier work 
show the behavior of relaxation times t q (obtained from 
NVE molecular dynamics of bulk equilibrium systems) as 
a function of £pps in all three systems. Reasonable data 
collapse is found, signifying a strong correlation between 
these quantities in all three systems. Fig. [3jd) shows a 
scaling plot where the independent variable takes an acti- 
vated form ln(r Q ) ~ £,pts/T, normalized by the values of 
the lengths and temperatures found from Fig.[3jb). The 
degree of collapse is similar to that shown in Fig. [3][c) 
with no adjustable parameters, making it difficult to dis- 
tinguish between scaling forms where the temperature 
dependence of the growth of r a may be attributed purely 
to the growth of £pts alone or the activated form com- 
monly found in the literature [H, 0, 13] ■ Furthermore, both 
Figs.[3](c) and[3](d) show slight but systematic deviations 
from perfect data collapse, suggesting that similar but 
distinct exponents associated with the scaling variable 
(£pts or ^pts/T) would be required to extend these 
scaling plots to larger t „ v alues; this may be similar to 
what is seen in Ref. [H [3]. Regardless, the correlation 
between £pts and relaxation time growth is striking. Fi- 
nally, we note that within the resolution of our current 
data, when relaxation times of these three models are 
similar, it is not only £pts a s defined here that matches, 
but rather the full overlap profile q(R); this is illustrated 
in the insets of Fig. O Hence, though our definition of 
£,pts is not unique, we expect any reasonable definition 
to give the same qualitative results as seen in Fig. [3) 
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FIG. 3. (a) Length vs. inverse temperature. For clarity, 
£ here represents £pts of the text. A clear distinction of 
scale between the WCA system and the other two systems 
should be noted, (b) The same data as in (a) but scaled to 
a temperature where all three models have the same length, 
T = 0.8 for the LJ and IPL and T = 0.6 for the WCA. 
Here the length £o is ~ 1.4. (c) Structural relaxation time 
T a vs. length £. The length scale correlates directly with 
the relaxation time scale, (d) Structural relaxation time vs. 
reduced length over reduced temperature. To and £o are the 
values giving collapse in (b). Error bars from our bootstrap 
analysis are approximately the size of the data symbols. All 
data are reported in the SI. 



The results presented here provide useful clues to the 
underlying causes of the viscous slowdown of dynam- 
ics during the supercooling process. £pts is the cen- 
tral structural feature associated with slow dynamics in 
RFOT-like theories [H, 0, 3 . In this sense, the results 
presented in Fig. [3] would seem to be in harmony with 
that viewpoint. It should be noted, however, that within 
the standard RFOT, the lengths in the relatively high 
temperature regime are not expected to necessarily cor- 
relate strongly with relaxation times, as they do here. 
Our results bear a resemblance to some features found 
in some short-ranged p-spin models, which serve as a 
possible paradigm for how the mean-field models that 
RFOT is based on are altered by fluctuations. In partic- 
ular, a weakly growing static length scale that is linearly 
correlated with Tln(r) all the way up through the high 
temperature regime as in Fig. |3]^d) are features found in 
the model of Ref. [H. Unfortunately, the static correla- 
tions in this model are spin-glass correlations, which are 
absent in real liquids. In general, studies such as those 
carried out here should shed light on how fluctuations 
modify (or destroy) the mean field behavior upon which 
the standard version of RFOT is based 2(| 27| . 

The mere existence of a growing £pts is not in contra- 
diction with a picture based on kinetically constrained 



models; however there are several aspects of our re- 
sults which seem difficult to reconcile with the facilitated 
picture. To address this, one has to examine results 
from thermodynamically interacting models which can 
be mapped onto kinetically constrained models. It has 
already been noted that local dynamics inside compact 
cavities show markedly increasing relaxation times as the 
cavity radius is decreased [HI]. While the same behav- 
ior holds in the three systems studied here (not shown), 
the opposite is true for large cavities in the triangular 
plaquette model, the one model known to be dual to a 
kinetically constrained model with hierarchical behavior 
0] ■ In addition, the absolute length scales found here are 
quite small. Since £pts might loosely be interpreted as 
the mean distance between defects in a facilitated model, 
this fact would seemingly translate into an unrealistically 
high density of defects with a rather weak temperature 
dependence. Indeed, the magnitude of £pts found in 
this work seems more consistent with the spatial extent 
of the "defects" found in a recent study than the separa- 
tion between these regions [28j |. Lastly, our results seem 
consistent with the idea that £pts an d £4 are distinct in 
the regime of weak to moderate supercooling |9(, which 
would be unexpected in models such as the triangular 
plaquette model [1, 0] It is important to bear in 

mind that plaquette models provide merely a small set 
of all possible mappings of an interacting system to one 
of kinetically constrained but thermodynamically non- 
interacting defects. It is unlikely, however, that any such 
mapping would have much to say about the small length 
scales found in this work, which are expected to be simi- 
lar to the coarse graining length scale in any mapping to 
a kinetically constrained model [3(| • 

In conclusion, £pts in three model glass forming sys- 
tems has been been extracted. Care has been exercised to 
avoid dynamical contamination of the overlap function. 
We find that £pts is small and grows systematically al- 
beit modestly as temperature is lowered. The spatial 
range of the thermodynamic overlap functions and ex- 
tracted £pts are clear discriminators of dynamical be- 
havior, even when simple structural features such as pair 
distribution functions are blind to differences in relax- 
ation times. We find that in the cases studied here, £pps 
correlates reasonably well in absolute terms with the re- 
laxation times in all three systems. Future work will be 
devoted to assessing if these correlations robustly hold 
over a wider range of temperatures and systems, as well 
as connecting the results uncovered here to other recent 
work on non-trivial structural and thermodynamic mark- 
ers of glassy behavior 0, [l2, 23 1. 
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SUPPLEMENTAL INFORMATION 



System Details 



chose to use box sizes approximately 1.2 times smaller, or 
/ 3 = 0.05. The bulk overlap for our systems is therefore 
pi 3 = 0.06. 



The LJ model is an 80:20 binary mixture of Lennard- 
Jones spheres with size and interaction parameters given 
by cr AB /<j AA = 0.8, (JbbIoaa = 0.88, and tAB^AA = 
1.5, cbb/caa — 0.5 [l5|]. The WCA system is charac- 
terized in the same way, but the interaction potential is 
shifted up by e y and cut off at r = 2 1 / 6 er ij - [16]. The IPL 
system is characterized by an inverse-power law poten- 
tial given by V{rij) = Aeij(a \j / 'r \j) n with n = 15.48 and 
A = 1.945 and the same size and energy parameters as 
the LJ system. All quantities are reported in standard 
reduced Lennard- Jones units, and all systems are studied 
at number density p = 1.2 [1 81 ] . 

Radial distribution functions have been generated for 
the three models at the temperatures studied in the main 
text. A figure showing these correlations for most of the 
temperatures in the relevant temperature regime can be 
found in Fig. |H Though differences do exist between the 
<?(r)'s at all temperatures, the differences stay relatively 
constant and are small. We do not attempt to demon- 
strate that these differences are too small to account for 
dynamical differences observed between the models, as 
this has already been discussed in prior work 17]. 



Simulation Details 

Eight independent configurations of 4050 particles 
were equilibrated at each temperature of interest for each 
system using Swap-MC [lo}. Six locations were chosen 
randomly in each of the eight independent configurations 
and particles whose centers fell outside a sphere of a given 
radius were frozen. This was done with the condition that 
the ratio of particle types and the density within a sphere 
of a given size matched that of the bulk. In order to en- 
sure this, radii were defined by the total number of par- 
ticles inside that spherical cavity i.e. R = (f^-) 1 ' 3 with 
N the number of particles and p the bulk density. Moves 
that took a particle outside of the cavity were rejected 
(although without this restriction, particles escaped only 
rarely and usually returned immediately). The overlap 
(averaged over the 48 independent cavities) was moni- 
tored as a function of time until reaching a plateau. This 
plateau value was extracted through a fit to a shifted 
stretched exponential and the long time limit was taken 
as the thermodynamic overlap. Error in this value was 
estimated by bootstrapping a distribution of curves to fit 

As in Ref. [13, we use 125 boxes (a 5x5x5 cube) to 
probe the overlap. In Ref. [l(J the authors choose / 3 = 
0.062876 which they found was small enough to prevent 
two particles from appearing in the same box at number 
density p = 1. Our systems studied were at p — 1.2 so we 



PSA Details 

In the Particle Size Annealing (PSA) method, the in- 
teraction between particles is adjusted such that er^- — >• 



an. Each fixed particle outside of 
= 1 and each cavity parti- 



2 ~~ 2 w y 
the cavity is assigned A 

cle is assigned a A^ based on some tunable switching 
function. The switching function we employ is given by 
X(t) = A + (1 - A ) [1 - (1 - (i/r,)) 7 ] where A is the 
smallest particle size, r s is the time over which switching 
occurs and 7 controls the steepness of the curve. For our 
purposes, we empirically chose to run at Ao = 0.6 for sev- 
eral bulk r a , and then to anneal with 7 = 4 and t s ten 
times longer than the length of the randomization phase. 
These simulations are run using the Swap-MC algorithm 
because for \{t) < 1, many swap moves are accepted 
and this greatly enhances sampling. Because these data 
are noisier than the thermodynamically averaged over- 
lap from an equilibrium simulation (due to the fact that 
only the 48 initial cavity configurations exist for refer- 
ence), two PSA simulations were run per cavity. These 
data were logarithmically binned and averaged across all 
of the cavities for a given system and temperature. The 
cavity overlap was taken by fitting the plateau of this 
curve to a horizontal line. 



Fits and Extraction of a Length 

Here we comment on the fitting form chosen for ex- 
tracting a length from the overlap data and report the re- 
sulting fit parameters. The authors of Ref. 10 find that a 

compressed exponential form q(R) = Aexp(— 
fits their low temperature data well. However, they fix 
7] = 1 to get good fits to their higher temperature data, 
resulting in a sharp jump in the lengths extracted. Phys- 
ically, we do not expect their fitting form (rather, any 
fitting form) to hold for all R > given that the over- 
lap quantity is not well defined for very small cavities. 
We expect that when R = 1, i.e. the cavity radius is 
approximately the size of a large particle " diameter" , to 
have a cavity with about one particle on average. For a 
cavity of this size, we would expect an overlap between 
zero and one, and moreover we predict this value to be 
approximately temperature independent and devoid of 
information related to emergent amorphous order. By 
still fitting our data to a compressed exponential, but 
shifting it to start at R = 1, it was found that the form 

q(R) = v4exp(- ) with a = 1 fits all of the data 

at all temperatures given a fixed A value in the range 
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FIG. 4. Radial distribution function for the three models. Shown are temperatures T = {0.9,0.8,0.7,0.6,0.55} for all models 
and also T = {0.5, 0.4} for the WCA model. Lines are shifted for clarity by 0.5 in panel (a) and by 0.75 in panels (b) and (c). 
The highest temperatures are at the bottom. 



(0.45 — 0.6). Thus a two parameter fit to all of the data 
becomes viable. By sweeping the value of A in this range, 
it was found that A = 0.5 fits the data best across the 
three models, though slightly better fits could be done 
by picking a different A for each model. We chose to 
use the same A value for all three models, because we 
felt this resulted in a value of £ that was more equivalent 
when comparing the models and because this decreases 
the possibility of spurious overfitting. From an alterna- 
tive perspective one sees by plugging in g(£ + 1) = 0.5/e 
that our fit form is reporting as a length the cavity size 
where q(R + 1) ~ 0.18 (again motivated by the physi- 
cal argument that nontrivia! amorphous order sets in at 
scales R > 1). Hence this definition could be used in- 
dependently of fit form if an alternative function were 
determined to fit the data better. 

Given this flexibility in method for extracting the 
length, we do not attempt to ascribe meaning to the fit 
parameter rj or to the exact magnitude of £pts- Never- 
theless, we show in Fig. [5] that generally the same quali- 
tative conclusions can be drawn for rj as for £pts, namely 
that the rj values can distinguish the LJ and IPL from the 
WCA at a given temperature but collapse when scaled 
to the onset temperatures. However, there is much more 
noise in the 77 values as extracted. As stated earlier, the 
errors in the fit parameters were determined by generat- 
ing a bootstrap distribution of q(R) curves from the 48 
standard q(R,t) or 96 PSA q(R,t) data and calculating 
the standard deviation of that distribution 



discussed in this section in Tab. Q] 



Computational Details 

The computations above comprised approximately 
30000 independent tasks broken into hundreds of thou- 
sands of jobs, totaling more than two million hours of 
computational time. They were executed in parallel on 

resources provided by the University of Chicago Compu- 
3.2 




32|. 



We report the fit values obtained using the methods 



FIG. 5. Fit parameter r] as a function of inverse temperature. 
Inset: 77 plotted against scaled inverse temperature where the 
scaling is as in Fig. 3(b). 



tation Institute, the Open Science Grid, and the Extreme 
Science and Engineering Discovery Environment. These 
resources were leveraged by expressing the simulations in 
the Swift parallel scripting language [31]. Example Swift 
code for this project can be found in Ref. [U 
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LJ 




IPL 




WCA 




Temp 




V 




V 




V 


1.0 


1.067 ±0.005 


1.599 ±0.019 


1.177 ±0.003 


1.667 ±0.012 










0.9 


1.243 ± 0.002 


1.853 ± 0.009 


1.305 ± 0.002 


1.787 ±0.010 


0.986 ± 


003 


1 


469 ±0.010 


0.8 


1.427 ±0.002 


2.008 ± 0.014 


1.466 ±0.001 


2.002 ± 0.000 


1.077 ± 


002 


1 


530 ± 0.006 


0.7 


1.620 ±0.002 


2.331 ±0.014 


1.581 ±0.119 


2.668 ± 0.984 


1.255 ±0 


003 


1 


778 ±0.011 


0.65 


1.770 ±0.003 


2.050 ±0.013 


1.865 ±0.002 


2.303 ±0.013 


1.322 ± 


002 


1 


783 ±0.007 


0.6 


1.949 ± 0.002 


2.329 ±0.011 


1.940 ± 0.003 


2.594 ±0.014 


1.423 ± 


003 


1 


858 ±0.011 


0.55 


2.055 ± 0.002 


3.119 ±0.017 


2.108 ±0.002 


3.332 ±0.016 


1.502 ± 


003 


2 


070 ±0.015 


0.5 










1.679 ± 


002 


1 


992 ± 0.001 


0.45 










1.884 ± 


002 


2 


755 ± 0.014 


0.4 










2.161 ±0 


003 


2 


527 ±0.013 



TABLE I. Values extracted from fitting as a function of temperature q(R) = 0.5 exp (^■^j^J f° r the three systems studied. 
Errors are calculated as discussed in the text of the supplemental information. 



